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Abstract. In fluorescence diffuse optical tomography (fDOT), the reconstruc¬ 
tion of the fluorophore concentration inside the target body is usually carried 
out using a normalized Born approximation model where the measured fluores¬ 
cent emission data is scaled by measured excitation data. One of the benefits 
of the model is that it can tolerate inaccuracy in the absorption and scattering 
distributions that are used in the construction of the forward model to some 
extent. In this paper, we employ the recently proposed Bayesian approxima¬ 
tion error approach to fDOT for compensating for the modeling errors caused 
by the inaccurately known optical properties of the target in combination with 
the normalized Born approximation model. The approach is evaluated using 
a simulated test case with different amount of error in the optical properties. 
The results show that the Bayesian approximation error approach improves 
the tolerance of fDOT imaging against modeling errors caused by inaccurately 
known absorption and scattering of the target. 


1. Introduction. Fluorescence diffuse optical tomography (fDOT) is an emerging 
imaging technique aiming at recovering the distribution of fluorophore marker inside 
diffusive target medium from measurements of fluorescent emission at the surface 
of the body [1, 2]. Typically, fluorescence agents bound to molecules or proteins are 
introduced to the bloodstream. The molecules then act as ligands as they attach 
themselves to the targeted receptor sites. The surface of the body is illuminated 
with light at the excitation wavelength of the fluorescence agent. The measurement 
system collects fluence from the body surface at the emission wavelength. The 
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inverse problem associated with fDOT is the estimation of spatially distributed 
fluorophore marker concentration in the body. 

In small animal studies [3, 2, 4], fDOT has been used for monitoring tumors in 
brain [5, 6], breast [7, 8] and lungs [9]. fDOT has also been used to monitor brain 
strokes [10], localizing lymph nodes that get affected at the onset of cancer [11] and 
studying the effects of drugs on tumor [12], In humans, fDOT has been used in 
imaging breast cancer [1, 13, 14, 15]. 

Computationally, the inverse problem in fDOT amounts to estimating the spa¬ 
tially distributed fluorophore concentration from the model 

V = A(na,n' a )h + e, (1) 


where h £ R" is the vector of unknowns, representing a pixel or voxel or other 
parametrization of the fluorophore concentration, y £ R m is the measurement vec¬ 
tor, e £ R m models the measurement noise and A(/z a ,/4) is matrix implementing 
the forward model (i.e., the light propagation model) corresponding to nominal ab¬ 
sorption and scattering distributions /i a and y ! s . The estimation of h from the model 
(1) is an ill-posed inverse problem, that is sensitive with respect to measurement 
and modeling errors. 

The fDOT inverse problem is most often carried out using the so-called normal¬ 
ized Born approximation model [16], where the measurement vector is the measured 
fluorescent emission data vector y £ bs scaled by the measured excitation data y® bs as 


V = 


yjbs 

2/obs’ 


( 2 ) 


and forward matrix A is 


-4(>a,/4) = diag(—^)i(> a ,/4) 

2/calc 


( 3 ) 


where A is the forward matrix corresponding to the raw fluorescence measurement 
and y® alc is computed excitation data corresponding to the absorption and scatter¬ 
ing distributions /z a and y! s [16, 17]. The convenience of the Born normalization 
comes from the fact that it does not require a reference excitation measurement 
from homogeneous reference media. The normalization also effectively calibrates 
the problem with respect to source strength and individual gains and coupling coef¬ 
ficients of individual source detector pairs [16, 4]. From the practical point of view, 
a further significant feature of the Born normalized model is that it can tolerate 
inaccurately known target absorption and scattering distributions (/x a , fi' s ) to some 
extent. 

The absorption and scattering distributions (/i a ,/4) are often interpreted as 
known (nuisance) parameters in the fDOT problem. However, in practical experi¬ 
ments, the actual values of these parameters are usually not known accurately, and 
therefore one is bound to use approximate measurement model 


y ~ A(n a ,*,li' s ,*)h + e, (4) 

where /r aj * and y,' s * are our estimates for the absorption and scattering of the body. 
Obviously, the use of incorrect realizations {/la.,*, a 4 *) in the measurement model 
induces unknown modeling error 

ai fJ’s) 

in the model, and since the inverse problem is ill-posed, this error may cause large 
artifacts to the reconstructed fluorophore image. 
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In this paper, we propose the reduction of the reconstruction errors caused by 
inaccurately known absorption and scattering of the body by the Bayesian ap¬ 
proximation error approach [18, 19]. We consider the approach using the Born 
normalized model (1-3), so that the starting point would be the model that has 
the best available tolerance for the inaccurately known absorption and scattering. 
The key idea in the Bayesian approximation error model is to represent not just the 
measurement error, but also computational model inaccuracy as a random variable. 
Hence, instead of the approximate measurement model (4), we write an accurate 
measurement model, 

y = +e (5) 

. - ' 

= + P (6) 

where the term e = [H(/i a ,/4) — H(/Lt a ,*,/4 *)]/i represents the approximation error 
and v = e + e the total noise. Obviously, the realization of e is unknown since it 
depends on the unknown h and the unknown nuisance parameters (/x a , /i(). However, 
in the Bayesian inversion paradigm, we can compute an estimate for the statistics of 
the approximation error e over the prior probability distributions of the unknowns 
and the nuisance parameters /z a and The approximation error statistics are then 
used in the inverse problem to compensate for the inaccurately known /r a and /i' s . 

The Bayesian approximation error approach was originally applied for discretiza¬ 
tion error in several different applications in Ref. [18]. For this reason, the term 
“approximation error” is commonly used also where “modeling error” might be a 
more appropriate term. The approach was verified with real EIT data in Ref. [20] , 
where the approach was employed for the compensation of discretization errors and 
the errors caused by inaccurately known height of the air-liquid surface in an indus¬ 
trial mixing tank. The application of the Bayesian approximation error approach 
for the discretization errors and the truncation of the computational domain was 
studied in Ref. [21], and for the linearization error in Ref. [22]. In Ref. [23] the 
approach was evaluated for the compensation of errors caused by coarse discretiza¬ 
tion, domain truncation and unknown contact impedances with real EIT data. In 
addition to EIT, the Bayesian approximation error approach has also been applied 
to other inverse problems and other types of (modeling) errors: Model reduction, 
domain truncation and unknown anisotropy structures in optical diffusion tomog¬ 
raphy were treated in Refs. [24], [25], [26] and [27]. Missing boundary data in the 
case of image processing was considered in Ref. [28]. In Ref. [29], again related to 
optical tomography, an approximative physical model (diffusion model instead of 
the radiative transfer model) was used for the forward problem. In Ref. [30], an 
unknown nuisance distributed parameter (scattering coefficient) was treated with 
the Bayesian approximation error approach. The compensation of errors caused by 
unknown optode coupling coefficients and locations was considered in Ref. [31]. The 
compensation of errors caused by unknown domain shape and discretization error 
was considered in Ref. [32]. The extension and application of the modeling error 
approach to time-dependent inverse problems was considered in Refs. [33], [34] and 
[35]. 

In this work, we evaluate the Bayesian approximation error approach with sim¬ 
ulated two dimensional (2D) and three dimensional (3D) examples. In the 2D 
example, we use different levels of severity of error in the nominal absorption and 
scattering distributions and show that the approach improves the tolerance of the 
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fDOT problem against inaccurately known absorption and scattering of the body 
over the conventional reconstruction approach. In the 3D example, we test the 
approach using simulated data using the Digamous atlas geometry [36, 37]. 

The remainder of the paper is organized as follows. In Section 2, the Bayesian 
approximation error model for the compensation of modeling errors due to unknown 
absorption and scattering in fDOT is presented. In Section 3, we review the light 
transport model which we use for the forward model. The details of data simulation, 
meshing, constructing the prior models and estimation of the approximation error 
statistics are described in Section 4. The simulation results are presented in Section 
5. Finally, the conclusions are given in Section 6. 


2. Bayesian approximation error approach. 

2.1. Statistical inversion in general. In the Bayesian approach to inverse prob¬ 
lems, the principle is that all unknowns and measured quantities are considered as 
random variables and the uncertainty of their values are encoded into their proba¬ 
bility distribution models [18, 38]. The complete model of the inverse problem is the 
posterior density model, that is, the information and uncertainities in the unknown 
or interesting variables given the measurements, given by the Bayes’ theorem 


7r(/i,/z a ,/4,e|y) 


ir(y\h, /z a , /z', e)n(h, p a , /z', e) 


n(y) 


( 7 ) 


where ir(y\h, /z a , p' s , e ) is the likelihood density modeling the probability of different 
measurement realizations when the realizations of h , /z a , a 4 an d e are given. The 
density n(h, /i a , p' s , e) is the prior model and it models our information on the un¬ 
known parameters before the actual measurements. The posterior (7) is practically 
always marginalized with respect to the unknown but uninteresting measurement 
related errors e as 

7r(ft,Ma,/4l2/) = J 7r ( /l .Ma,Ms,e|?/)de, (8) 

for details in the case of the additive error model y = A(h) + e, see Refs. [30] and 
[39], for other types of errors see Ref. [18]. The posterior density n(h, /z a , y! s \y) is a 
probability density in a very high-dimensional space. Thus, in order to get practical 
estimates for the unknowns and visualize the solution, one needs to compute point 
estimate(s) from the posterior density, the most common choice in high dimensional 
cases being the maximum a posteriori (MAP) estimate. In principle, one could 
attempt to compute the MAP estimate for all the unknown model parameters 


( h ,/z a , /z')map = arg max tt(/i,/ z a ,/z'| y). (9) 

See Refs. [40], [41] and [42] for simultaneous reconstruction of fluorescent and optical 
parameters. Also, in Ref. [43] a method to determine a modified optical parameter 
(related to /z a , /z') from measured excitation data and include it as apriori anatom¬ 
ical information into fDOT imaging is presented. Alternatively, one could treat 
the uncertainty in the values of nuisance parameters (/z a ,/z') by marginalizing the 
posterior density as 


tt(%) 




( 10 ) 


and then compute estimate for the primary unknowns from the posterior ir{h\y). 
However, the solution of (10) would require Markov chain Monte Carlo integrations 
that would be computationally infeasible for practical purposes. 
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The key idea in the Bayesian approximation error approach is to find an ap¬ 
proximation n(h\y) for the posterior (10) such that the marginalization over the 
uncertainty in the values of (/z a , n' s ) is carried out approximately and in a computa¬ 
tionally feasible way. 

Before presenting the Bayesian approximation error approach for treating the 
uncertainty in the optical parameters (p, a ,t4)> we first review the standard fDOT 
reconstruction approach where (/x a , p' s ) = (pL a ^,p' s *) are treated as known and fixed 
variables. 


2.2. Conventional error model. In most of fDOT reconstruction schemes, the 
optical parameters (/r a ,/4) are treated as known (fixed) nuisance parameters with 
values p a = p aj * and p' s = p' s *. Within the Bayesian setup, this is equivalent to 
considering (/x a , n' s ) as fixed conditioning parameters, leading to posterior density 
model 


Tr(h,e\y,p a = /z a ,*,/4 = //',*) 


*{y\h, = // a ,*, /4 = Ms,*) e)7r(fe)7r(e) 

Ay) 


( 11 ) 


However, in cases that the values (/z ai *,/4*) are incorrect, the model (11) can be 
grossly misleading. 

Given the observation model (4) with fixed realizations (/r a ,/4) = (p a ,*,/4,*)i 
and modeling the measurement noise as normal e ~ A/"(e»,T e ), where e* E K. m is 
the measurement noise mean and T e E M mx " 1 is the measurement noise covariance 
and marginalizing (11) over the unknown measurement errors e as 


a,* 



Ah, e|y,/r a ,*, 


/4,*) de 


the posterior density becomes [30, 39] 


7r(%,Ma,*,/4,*) « exp |-^l|y - 


- e* 


r (h) 


( 12 ) 


In addition, if the random measurement noise has zero mean (e» = 0), the MAP 
estimate corresponding to the posterior (12) is obtained as 


hcem = arg max n(h\y, p a = /x a ,*, /4 = p' s *) 

h 

= arg min{||i/ —A(/i a *,/i' )/i||p_i — 21og7r(/i)}, (13) 

h ’ 1 e 


where the Cholesky factor LjL e = P e x . We refer to the solution of (13) as the 
MAP estimate with the conventional error model (CEM) approach. 


2.3. Bayesian approximation error model. In the Bayesian approximation er¬ 
ror approach, we write the measurement model as 

V = A(n a ^,fj,' s ^)h+ [A(p, a , p,' s )h - A(p, a ^,p! s ^)h] +e (14) 

-v-' 

eOaX.M 

= A^^^sJh + v (15) 

Note that the model (15) is exact, see equation (1). Here the term e = [A(/r a , fi' a ) — 
A(/x a> *, a4 *)]/i represents the approximation error and v = e + e the total error. 
Obviously, the realization of e is unknown since it depends on the unknowns h 
and the unknown nuisance parameters (p a ,p,' s ). However, in the Bayesian inversion 
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paradigm we can compute an estimate of the statistics of e using the prior probabil¬ 
ity distributions of the unknowns and the nuisance parameters and approximately 
marginalize over e. 

To include the uncertainty in the noise process e into our computational posterior 
model, the core step in the Bayesian approximation error approach is approximate 
pre-marginalization of the joint distribution of the parameters ( y , h, e, e, /r a , /i' s ) over 
the nuisance parameters (e, e, p, a , /z'). Following the approach in Refs. [30] and [39] 
and making a Gaussian approximation for the joint density tt( h,e), we obtain the 


approximate likelihood model 



Tr(y\h) oc exp | - 

\h - ,*) h - V*\hWl-i h | 

(16) 

where 



^*\h 

e * + £ * + r e,hF h 1 (ft - K) 

(17) 

Ey| h = 

r e + r E - 

(18) 


For the posterior model, we write an approximation n(h\y) oc Tr{y\h)n(h) 1 leading 
to MAP estimation problem 


Ke m = argmin{||y - A(/z a *,/4 *)/i - v*\ h ||p,_i -2 log7r(/i)}. (19) 

We refer to the estimate (19) as the MAP with the Bayesian approximation error 
model (MAP-AEM). 

In the following, we employ the enhanced error model [18, 19], where we approx¬ 
imate the error e and the unknown h as mutually uncorrelated, implying that the 
mean and covariance in equations (17)-(18) become 

^*i?i ~ e* + £*, v v \ h « r e + r e . 

In this work, the mean e* and covariance r e were computed by sampling based 
Monte Carlo integration (see Section 4.3). 

2.4. Estimates to be computed. To evaluate the Bayesian approximation error 
approach, we compute the following estimates 

(MAP-REF): The unknown h ce m using correct fixed optical coefficients (/x a , /4) : 

h le f = argmin{||y - A(/z a , fJ.' a )h\\i-i - 2log tt(/i)}. (20) 

In other words, in MAP-REF (/x a ,/x') are known exactly. Therefore, this 
estimate serves as reference of conventional reconstruction when (/x a , /i' s ) are 
known. 

(MAP-CEM): The unknown h cem using incorrect fixed optical coefficients 
(/Xa,*,A*s*) (Eq. 13). This estimate serves as a conventional reconstruction 
when nominal (/x a ,/x') are incorrect. 

(MAP-AEM): The unknown h aem using the same fixed incorrect coefficients 
(/Xa,*,/4 *) (Eq. 19). This estimate serves as a Bayesian approximation error 
model reconstruction when nominal (/r a ,/x') are incorrect. 

The values (/z a ,*, *) i n the computations correspond to the expectations of the 

(proper Gaussian smoothness) prior models 7r(/x a ) and ), respectively. The prior 
means were homogeneous distributions with y a ^(r) = O.Olmm -1 and y' s *(r) = 
lmm” 1 . 
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3. Forward model. Let D C R n , n = 2,3, denote the object domain. In a diffusive 
medium like soft tissue, the commonly used light transport model for excitation and 
fluorescence light is the diffusion approximation (DA) to the radiative transport 
equation (RTE) [44]. In this paper the DC (zero-frequency) domain version of the 
diffusion approximation is used 

(—V ■ n{r)V + /r a (r)) $ e (r) = 0, r G D, (21) 


, e / \ 1 , . <9<h e (r) 

$ (r) + 



r G r s 
r G \ r s 


( 22 ) 


(—V • «(r)V + Hair)) $ f (r) = /i(r)4> e (r), reH, (23) 


$f ( r ) + = °’ ( 24 ) 

where <f> e (r) := <h e is the excitation photon density, 4> f (r) := 4> f is the fluorescence 
emission photon density, /z a (r) := /z a is the absorption coefficient, re(r) := k is the 
diffusion coefficient. The diffusion coefficient k is given by k (r) = 1 /(n(/i a (r) + 
A t s( r )))> where /4(r) := /x' is the reduced scattering coefficient. For simplicity 
the spectral dependency of the optical properties (/z a ,/z') is omitted and they are 
modelled the same at the excitation and emission wavelengths. The parameter 
h(r) := h is the fluorophore concentration. The parameter q(r) is the strength of 
the light source at location r s C dQ. The parameter £ is a dimension dependent 
constant (£ = 1/ir when JlcK 2 ,( = l/2 when ficK 3 ), aisa parameter governing 
the internal reflection at the boundary dQ and i) is the outward normal to the 
boundary at point r. The measurable excitation data and fluorescence emission 
data (Eq. 2) are given by 


V e {r) 


V l {r ) 


a$ e (r) 2y 

= — &(r), 

ov a 

9$ f (r) 2 7 f 

-« = —$ (r), 

cw a 


t G r d , 


r G r d , 


(25) 

(26) 


where C are the detector locations. 

For the inverse problem, the mapping A(/x a ,/z')/i (Eq. (1)) is given by [16, 4, 17] 


A(// a , n' s )h 


f n ® e (r s ,r)V e (r d ,r)h(r)dr 

fn 4> e (r s ,r)dr 


(27) 


where $ e (r s ,r) is the computed excitation photon density due to source q(r). 
d/ e (r d ,r) is the computed adjoint solution (photon density due to sources placed at 
detector locations r d ). 

The numerical approximation of the forward model used here is based on a finite 
element method (FEM) solution of Eq. (21-26). 


4. Computation details. 


4.1. Simulation of the measurement data. 
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4.1.1. 2D simulations. In the 2D numerical studies, the domain SlcM 2 was a disk 
with radius r = 25mm. The measurement setup consisted of 16 sources and 16 
detectors. The source and detector optodes were modeled as 1mm wide surface 
patches located at equi-spaced angular intervals on the boundary dD. With this 
setup, the vector of fDOT measurements (2) was y G R 256 . Five targets with 
different optical properties of absorption and scattering were used in the simulations 
(see Fig. 2, column 1, 2). The fluorophore concentrations of these targets is shown 
in column 1, Fig. 3. For the simulation of the measurement data, a mesh with 
33806 nodes and 67098 triangular elements was used. In the inverse problem, a 
FEM mesh with 26075 nodes and 51636 elements was used. 

4.1.2. 3D simulations. In the 3D simulations, the domain Q C R 3 was a mouse 
atlas “Digimouse” [36, 37]. The measurement setup consisted of 32 sources and 32 
detectors (see Fig. 4, bottom right). The source and detector optodes were modeled 
as 1mm wide surface patches placed on the top-surface of the boundary dQ. With 
this setup, the vector of fDOT measurements (2) was y £ R 512 . The optical prop¬ 
erties of the target are shown in column 1, Fig. 4. The fluorophore concentrations 
of the target is shown in column 2, Fig. 4. For generating the measurement data 
and for the inverse problem we used the same mesh obtained from the Digimouse 
website [37] which had 58244 nodes and 306773 elements. 

The simulated measurement data was generated using the FEM approximation 
of Eq (21)-(26). Random measurement noise, drawn from a zero-mean Gaussian 
distribution was added separately to the simulated measurement excitation and 
fluorescence emission data as, 



(28) 

(29) 


Here e f ~ A/"(0, diag(cr e f >1 ,.., cr e t m )) £ R m is the noise added to the fluorescence 
emission data, e e ~ jV(0, diag(cr e e , 1 , • cr e e ,m)) £ R m is the noise added to the 
excitation data. The standard deviations {a e e > i, ..,<T e e m } and {cr e f ,i> ••i <T e f ,m} were 
specified as 1% of the simulated noise free measurement excitation and fluorescence 
emission data, y® alc and y l calc respectively. 

In order to estimate the noise statistics for the inverse problem, a Gaussian 
approximation for the measurement noise e, 7r(e) = A/"(0,r e ) was constructed. The 
covariance r e was approximated by a diagonal where the diagonal elements (i.e., 
standard deviations of the measurements) were estimated as sample averages from 
100 noisy realizations of the Born ratio (Eq. (28)). In a practical setup, similar 
estimation of T e can be carried out by repeated measurements from a phantom 
target. 

4.2. Prior models. For drawing the samples of optical coefficients x, where 


x = (/T a ,/x',/i) T , 


for the estimation of the approximation error statistics, we used a Gaussian (Markov 
random held) [45] prior model along with non-negativity constraint [46], 


7 t(x) = ttg(x)tt + (x). 


( 30 ) 
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Here 'xg(x) is a Gaussian smoothness prior model and 7r + (x) is the non-negativity 
constraint. Section 4.2.1 describes the implementation of the Gaussian smooth¬ 
ness prior model, Section 4.2.2 describes the implementation of non-negativity con¬ 
straint. 

4.2.1. Proper Gaussian smoothness prior. As we need to draw samples of the un¬ 
known x for the estimation of approximation errors, we need a proper (integrable) 
prior distribution for the unknowns. In this study we used a proper Gaussian 
smoothness prior 710 ( 2 ;) for the unknowns. In this model, the absorption, scatter¬ 
ing and fluorophore concentration images p a , p' s and h were modeled as mutually 
independent Gaussian random fields with a joint prior model 

7r G (a;) oc exp{ —- £*)|| 2 }, L^L X =T~ 1 (31) 

where 



( Ma,* \ 

1 

/ IV 

0 

0 \ 

X* = 


, r* = 

0 

r* 

0 


V k y 

> 

V 0 

0 

V / 


In the construction of the mean vectors p a ,*,^ „,/i» and covariances T^, and 
Th, the random field, say / (i.e. , either p a , p' s or h), is considered in the form 

/ = /in + /bg 

where f ln is a spatially inhomogeneous parameter with zero mean, 

/ in ~ AT{0, T inJ ) 

and /b g is a spatially constant (background) parameter with non-zero mean. For 
the latter, we can write /b g = qf, where I £ 1" is a vector of ones and q is a 
scalar random variable with distribution q ~ A /"(c, <rj( g ^). In the construction of 
rin,/, the approximate correlation length can be adjusted to match the size of the 
expected inhomogeneities and the marginal variances of /&:s are tuned based on 
the expected contrast of the inclusions. We model the distributions /j n and /b g as 
mutually independent, that is, the background is mutually independent with the 
inhomogeneities. Thus, we have 

f* = cl, Tf = F in J + al gJ U T 

See [18, 24, 30] for further details, and see [47] for an alternative construction of a 
proper smoothness prior. 

The parameters in the prior model tt(x) were selected as follows. The mean 
for background absorption, scattering and fluorophore concentration were set as 
= 0.01mm -1 , p' s * = 1mm -1 and h* = 0mm -1 . The standard deviations 
<Tb g ,/i a and CTb glM ' of the background values were chosen such that 2 s.t.d. limits 
equaled 25% of the mean values p a ,* and p' s * and the standard deviation <7b g,h 
of the background value for h was chosen as 0.25 mm -1 . In the construction of 
Tin,/, the correlation length for p a , p' s and h in the prior was set as 16mm. The 
marginal standard deviations were set to equal values in each pixel and <Ti n , Ma and 
were chosen such that 2 s.t.d. limits equaled 50% of the mean values p a ,* 
and p' s *. a m ,h was chosen as 1mm -1 . Thus, the overall marginal variances (i.e., 
diagonal elements of T Ma , r M ' and T^) were 2cr Ma = 0.0056mm -1 , 2tr M / = 0.56mm -1 
and 2ah = 1mm -1 . This gives overall 2 s.t.d. intervals p a £ [0.0044,0.0156]mm -1 , 
p' s £ [0.44,1.56]mm -1 and h £ [—1,1] i.e., the values of absorption, scattering 
and fluorophore concentration are expected to lie within theses intervals with prior 
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probability of 95%. The same prior parameters (correlation length, means and 
standard deviations) are used in 2D and 3D priors. 

4.2.2. Non-negativity constraint. In addition to the Gaussian smoothness prior, a 
non-negativity constraint was applied while drawing samples for the computation 
of approximation error statistics. A tolerance value for the values of x was chosen 
and the values of x drawn from prior ttg{x) that were less than the tolerance value, 
were set equal to the tolerance value as, 

x = max(x,tol). (32) 

Here “tol” is the tolerance value, tol = 10e -6 was used in this study. 

Two random draws of x = (/z a , fi' a , h) T from tt(x) in the 2D simulation domain 
are shown in Figure 1. 




0 0.02 0 2 0 1 

Figure 1. Two draws (top and bottom row) of /x a , /z' and h from 
w(x) in the 2D simulation domain. From left: first column is /z a , 
second column is g! s and third column is h. 

The non-negativity prior during the computation of the MAP estimates (20) 
, (13), (19) is taken into account by applying an exterior point search method 
[46]. In this method, the non-negativity problem is approximated by a sequence of 
unconstrained problems 

hMAP = argminj —log(7r(%)) — log(7r(ft))}, (33) 

h 

= argmin{—21og(7r(%)) + \\L h {h - h*)\\ 2 + E J (h)} (34) 
h 

where Tr(h\y) is the likelihood and 7r(h) is the prior probability distribution. E° (h) 
is a penalty functional that is used to penalize the negative components of the 
solution h, the super index j is used to denote the j th problem in the sequence. 
The mean h * G 1" and Cholesky factor L^Lh = G R raxn is obtained from the 
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Gaussian smoothness prior model 7 rc(a;). In this study we employ a functional of 
the form 

n 

E j (h) = E 7 J (/>(hk) (35) 

k =1 


where 


<t>(h k ) 


( hk ) 2 , hk < 0 
0 , otherwise, 


(36) 


and 7 3 ,j = 1,2, ,.,M is a sequence of increasing positive penalty parameters. The 
exterior point methods guarantee the non-negativity of the solution only in the 
asymptotic limit j —>• oo. In this paper, we used { 7 1 = 1 , 7 2 = 10, 7 3 = 100} as 
the sequence of exterior point parameters. The incorporation of the non-negativity 
constraint to the MAP estimates (20) , (13), (19), leads to non-linear minimization 
problems, which were solved by a Gauss-Newton algorithm with an explicit line 
search algorithm [48]. 


4.3. Estimation of approximation error statistics. For the computation of 
the error statistics, first the samples of absorption, scattering and fluorophore con¬ 
centrations 


S = {x^\ l = 1,.., N s } 


(37) 


were drawn. Two of the samples are shown in Figure 1. The samples were then 
used for the computation of the accurate forward solutions A (^and ap¬ 
proximate forward solutions A(/z ai *, fi' s *)h^\ and the samples of the approximation 
error 


e^ = [A(^\^)-A(^,^)]h^\ 

were computed. Then the mean e* and covariance r e were estimated using the 
samples {e^} as 




1 

n7 


N s 

Ye® 


e =1 


(38) 


r e 



l ' 1 s 

E( ffW 


-e*)(e (f) -e *) 1 


(39) 


5. Results. 

5.1. 2D simulations. The true absorption and scattering parameters (// a , fi' s ) used 
in the simulations are shown in Figure 2, column 1 (absorption) and column 2 
(scattering). The nominal values (/i a ,*j/4*) that are used in the estimates h cem 
and /laem are shown in Figure 2, column 3 (absorption), column 4 (scattering). The 
true target fluorophore distribution h is shown in first column in Figure 3. The 
MAP estimates with the measurement data from the target domains in columns 1 
and 2 in Figure 2 are shown in Figure 3, column 2-4, in matching order of rows. 
The estimates are: 

(MAP-REF): The MAP-ref estimates using correct fixed (/x a , fj,' s ) 

href = argmin{||y - A(// a , /4)/i||p_ 1 + \\L h (h - /i*)|| 2 + E 3 (h)}, (40) 

h 1 e 

are shown in column 2, Figure 3. This corresponds to the reference estimate 
of conventional reconstruction when (/i a ,/r() are known exactly. 
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(MAP-CEM): The MAP-CEM estimates with fixed optical coefficients (/z aj *, /z' * 

hcem = argmin{||y - A(/z a ,*,/z' *)/i|| 2 i + \\L h (h - h *)|| 2 + E J (h)}, (41) 

are shown in column 3, Figure 3. This corresponds to estimate with conven¬ 
tional reconstruction when the nominal values of (/z a ,/z') are incorrect. 
(MAP-AEM): The MAP-AEM estimates with the same fixed optical coeffi¬ 
cients (/Za,*,^,*) 

h ae m = argmin{||?/ - A(/z a ,»,/z' *)/i - zz*|fc|| 2 i + \\L h (h - /i*)|| 2 + E J (h)}, (42) 

h ’ 1 u\h 

are shown in column 4, Figure 3. This corresponds to estimate with Bayesian 
approximation error model when the nominal values of (/z a ,/Zg) are incorrect. 
The relative error in the MAP estimates (40), (41), (42), 

Error = x 100%, (43) 

|| 'Hrue || 

where h is the estimated fluorophore distribution and ht rue is the true fluorophore 
distribution are shown in Table 1. 

Table 1. Relative error (%) in MAP estimates (40), (41), (42) for 
each 2D simulation test cases (test cases are numbered from top to 
bottom in Fig 1 and 2). 


Case 

^-ref 

^cem 

^aem 

1 

42 

42 

59 

2 

38 

64 

61 

3 

44 

100 

66 

4 

42 

117 

62 

5 

43 

116 

66 


5.2. 3D simulations. The true absorption and scattering parameters (/z a , /z') from 
horizontal and vertical slices of target mouse model used in the simulations are 
shown in Figure 4, column 1: top (absorption) and bottom (scattering). The true 
target fluorophore distribution h is shown in second column in Figure 4. The nomi¬ 
nal values (p, ai *, /z' *) that are used in the estimates h cem and h aem are homogeneous 
distributions with = 0.01mm -1 and /z' *(r) = 1mm -1 . The MAP estimates 

are: 

(MAP-REF): shown in column 3, Figure 4. This corresponds to the reference 
estimate of conventional reconstruction when (/z a ,/z') are known exactly. 

(MAP-CEM): shown in column 4, Figure 4. This corresponds to estimate with 
conventional reconstruction when the nominal values of (/z a ,/z') are incorrect. 

(MAP-AEM): shown in column 5, Figure 4. This corresponds to estimate 
with Bayesian approximation error model when the nominal values of (/z a , /i') 
are incorrect. 

The results show that the Bayesian approximation error model efficiently com¬ 
pensates for the modeling errors caused by inaccurately known (/z a , /z') of the body 
and produce estimates that are qualitatively similar to the conventional estimates 
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|T =0 

^ s 

LI =0 
a 



|LL’ =2 

^ s 

|i -0.02 


Figure 2. First and second columns show /z a (left) and /z' (right) 
of the body Q in test cases 1-5 (top to bottom). Third and fourth 
columns show the (incorrect) absorption and scattering /z a ,* and 
/i' * that are used the in the computation of the estimates MAP- 
CEM and MAP-AEM. 


with exactly known /z a and /z'. The inclusions in the fluorescent contrast are well lo¬ 
calized with the Bayesian approximation error model and the estimates are relatively 
free of the artifacts that are present in the estimates h cem with the conventional 
error model using the same incorrect values of (/z aj *,/z' *), see rows 2-5 in Figure 3 
and columns 3-5 in Figure 4. 
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h (true) MAP-ref MAP-CEM MAP-AEM 




0 1 


Figure 3. First column: fluorophore distribution h of the body 
ft, in test cases 1-5 (top to bottom). The second column shows 
the MAP-ref estimate using the correct absorption and scattering 
(forward matrix A(/i a ,/4))- The absorption and scattering images 
/t a and /t' are shown in columns 1 and 2 in Figure 2 (rows in re¬ 
spective order). Third column shows the the MAP-CEM estimate 
using the incorrect absorption and scattering values (forward ma¬ 
trix ^4(/x a ,*, yUg *)). The /r ai * and /i' * are shown in third and fourth 
column in Figure 2. The fourth column shows the MAP-AEM 
estimates using the same incorrect forward matrix A(/t a ,*,/t' *). 

A notable feature in the reconstructions is that estimate of h has lower contrast in 
the reconstructions with the approximation error model than in the reconstructions 
with the conventional noise model. This can be explained by the fact that covari¬ 
ance of the combined noise e + e is larger than the covariance of random noise (i.e., 
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0 2 

Figure 4. First column: target \x a (top) and /x' s (bottom) at one 
horizontal and vertical slices of 3D mouse model. The second col¬ 
umn shows the respective slices target fluorophore distribution h 
of the body. Third column shows MAP-ref estimate using the cor¬ 
rect absorption and scattering values (forward matrix A(/z a , /z()). 

Fourth column shows MAP-CEM estimate using the incorrect ab¬ 
sorption and scattering values (forward matrix A(/z a *, n' s *)). Fifth 
column shows MAP-AEM estimates using the same incorrect for¬ 
ward matrix A(/z a> *, fx' s *). Bottom-right: Mouse surface with po¬ 
sition of sources marked with yellow asterisk(*) and position of 
detectors marked with blue asterisks. 

r e + r e > P e ), implying that the relative weight of the data residual term compared 
to the prior (or regularization) term becomes smaller in the MAP estimate with 
the approximation error model compared to the conventional noise model. Thus, 
the estimate gets, loosely speaking, drawn more strongly towards the prior mean, 
leading to a loss of contrast. This is also seen from the error estimates in Table 
1. In the first row which corresponds to the case that (/z a ,*,l4*) are correct, the 
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estimation error with the approximation error model is larger than with the con¬ 
ventional noise, and this discrepancy in the error arises from the lower contrast in 
the estimate of h with the approximation error model. However, when the values 
(Ma,*)A4*) are incorrect, rows 2-5 in Table 1 and Figure 3, the estimation error 
with the approximation error model is smaller than with the conventional noise 
model, and moreover, the estimation error does not change much as the distance of 
g *) from the true values (/z a> A t s) increases when moving from row 2 to row 
5 in the Table and Figure. 

In Figure 4 we show simulation using a realistic 3D mouse model where we 
simulated high accumulation of fluorophores in the kidneys and lungs of the mouse. 
To make the simulation physiologically realistic we simulated background physiology 
by adding random fluorophore concentration h = ma x(h', 0) where h' ~ Af( 0, cr 2 I) 
with a = 0.2mm , all over the mouse domain except the lungs and kidneys. The 

localisation of fluorophores in the kidney and lung positions can be seen in the 
MAP-ref estimates (column 3 Figure 4). They appear slightly distorted in the MAP- 
CEM estimates (column 4 Figure 4). However, the fluorophore concentrations are 
relatively better localised in the MAP-AEM estimate (column 5 Figure 4) with a 
slight loss of contrast. 

6. Conclusions. In this paper the recovery from errors caused by incorrectly mod¬ 
eled absorption and scattering in fDOT was considered. Born ratio is known to toler¬ 
ate artefacts due to unknown absorption and scattering to some extent. However, in 
case the absorption and scattering properties are highly heterogeneous, incorrectly 
modelled absorption and scattering induces errors in the fDOT reconstructions. 

In this paper, the Bayesian approximation error approach was applied for the 
compensation of the errors caused by unknown absorption and scattering in fDOT. 
The modeling errors caused by the inaccurately known absorption and scattering 
were modeled as an additive modeling error noise in the observation model, and the 
posterior density model was then marginalized approximately over the unknown 
modeling errors by using a Gaussian approximation for the joint statistics of the 
primary unknown (fluorophore concentration) and the modeling errors. 

The approach was tested with 2D simulations with various target distributions of 
absorption and scattering. The results show that the approximation error model can 
efficiently compensate for the reconstruction artefacts caused by unknown absorp¬ 
tion and scattering coefficients, even in the cases of highly heterogeneous absorption 
and scattering coefficients where the conventional estimates using the Born ratio 
contained severe artefacts. The approach was also tested with a 3D simulation 
using a mouse atlas. The MAP estimates using the Bayesian approximation error 
model show better localisation of the fluorophore concentration compared to the 
conventional estimate with the normalised Born approximation model. 

A tradeoff of the Bayesian approximation error model was found to be a small 
loss in contrast of the estimated fluorophore concentrations. We suggest that the 
approximation error model would be a feasible complement to the Born ratio model 
for handling the uncertainty related to the unknown absorption and scattering pa¬ 
rameters in fDOT. 
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